This invention relates to the field of geophysical prospecting and, more particularly, to a method for determining refraction static solutions to apply to seismic data.
In the oil and gas industry, geophysical prospecting techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a seismic energy source is used to generate a seismic signal which propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflections are recorded by seismic detectors located at or near the surface of the earth, in a body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the location of the subsurface reflectors and the physical properties of the subsurface formations.
In the FIG. 1, a cross-sectional view of the earth that includes a weathering layer 10 beneath the earth""s surface E is shown in the vicinity of a subsurface seismic survey. The weathered layer 10 is often referred to the low-velocity-layer or LVL. Seismic energy is imparted into the earth at a shotpoint S, typically at the surface or some suitable depth within the weathering layer 10. The seismic energy travels outwardly from the shotpoint S through the weathering layer 10 and therefrom through an interface I between the weathering layer 10 and deeper subsurface earth formations 16. The seismic energy travels into the earth and portions of the energy are reflected by interfaces, such as one indicated at 14. Reflected energy travels upwardly from the reflectors through the subsurface formations 16 through the weathering layer 10 to a detector or array of detectors indicated schematically at R located at or near the earth""s surface E. The responses of the reflectors are then recorded and processed. It should also be understood that the present invention may be used with marine seismic survey data, as well. The purpose of reflection surveying is to identify subsurface formations or features of interest.
FIG. 1 contains a seismic datum labeled Datum. The seismic datum is an arbitrary reference surface, the reduction to which minimizes local topographic and near-surface effects. Seismic times and velocity determinations are referred to the datum plane (usually, but not necessarily, horizontal and planar) as if sources and geophones had been located on the datum plane and as if no weathered layer or low-velocity layer existed.
However, the low-velocity or weathering layer 10 is not of uniform thickness or elevation and the composition of the materials and their density within the layer varies, as does the seismic velocity within the weathering layer. The lack of uniformity in weathering layer characteristics introduces unwanted effects, anomalies and fluctuations in seismic energy traveltimes. These difference in seismic energy traveltimes can induce severe distortions in seismic images during data processing. These effects are generally referred to as statics, and it is well known in the art that static corrections may be applied to the data to correct for these effects.
FIG. 1 shows an example seismic raypath P1 through P4 traversing the weathered layer 10 and the subsurface earth formation 16. The time seismic energy takes to traverse the weathered layer may be highly variable both in traversing from the source S through the weathered layer 10 along P1 as well as along P2, P3 and returning to the detector R on the surface along P4.
Static corrections are corrections applied to seismic data to compensate for the effects of variations in elevation, near-surface low-velocity-layer weathering thickness, weathering velocity, and/or reference to a datum (See Sheriff, R. E., 2002, Encyclopedic Dictionary of Exploration Geophysics: Soc. of Expl. Geophys., 334-335). The objective is to determine the reflection arrival times that would have been observed if all seismic recording measurements had been made on a flat plane with no weathering or low-velocity material present, equivalent to the situation if the source S and receiver R were positioned at the level of the Datum of FIG. 1. After static corrections are applied, the source S would appear as source Sxe2x80x2 at the elevation of the Datum, and receiver R would appear as Rxe2x80x2 at the elevation of the seismic datum. These corrections are based on uphole data, refraction first-breaks, event smoothing, and sometimes other geophysical methods. The most common convention is that a negative static correction reduces the reflection time. Uphole-based statics involve the direct measurement of vertical traveltimes from a buried source. This is usually the best static-correction method where feasible.
Underlying the concept of conventional static corrections is the assumption that a simple time shift of an entire seismic trace will yield the seismic record that would have been observed if the geophones had been displaced vertically downward (or upward) to the reference Datum, an assumption not strictly true, especially if the surface-to-datum distance is large, and that the velocity of the weathering layer 10 does not change horizontally. Conventional static correction methods are most apt to fail where there are 1) large rapid changes in the topography or base of weathering, 2) horizontal velocity changes below the weathering, thus violating the assumption that the subdatum velocity does not vary significantly, 3) large elevation differences between the datum and the base of the weathering, or 4) inadequate controls on long-wavelength statics. Large sea-floor relief is apt to be associated with horizontal velocity changes that cannot be compensated with static corrections.
In seismic data processing, it is desirable to correct for these statics and eliminate as much as possible the effects of variations in the weathering layer 10 and other statics on the seismic data. It is highly desirable to have the seismic data be in a form as if it had resulted from a survey conducted on a substantially flat plane or datum at a constant elevation in the earth. To compensate for statics, it is necessary to determine the amount of time delay introduced by travel of seismic energy above or below the datum level and then remove the effects of this time delay from the seismic data. In effect, static correction compensates for elevation differences of the source S and detector R from the datum level, for changes in thickness in and along the weathering layer, and for variations in the density and velocity of the weathering layer 10.
If the structure and dynamics of weathering layer anomalies were known, the best way to tackle this problem would be to perform wave-equation datuming or depth migration from the surface through the known structure. However, 3-D prestack depth imaging and datuming are computationally expensive. Therefore, prior art statics applications that assume surface-consistent ray propagation through the near surface weathered layer have remained the main tool to account for near-surface anomalies.
Prior art refraction statics solution algorithms have generally employed a four step procedure. 1) 1) First arrival times of refracted arrivals are determined and properties of the weathered layer (such as weathered layer velocity) are estimated. 2) A trial model of the weathered layer is developed. This model may consist of one or more layers, or for tomographic algorithm approaches, a large number of small cells. 3) An iterative scheme is employed to adjust the model parameters so that they become more consistent, in some optimum sense, with the first arrival timing information and other user specified information. Many iterative schemes have been used to invert for model parameters, for example a generalized linear inversion (GLI) scheme or a tomographic scheme (Hampson, D. and Russell, B., 1984, First-break interpretation using generalized linear inversion, 54th Ann. Internat. Mtg: Soc. of Expl. Geophys., Session:S10.4). 4) Vertical rays are traced through the final version of the model at source and receiver locations to determine the static corrections that can be applied to the seismic data.
These traditional methods have several shortcomings. First, these methods require an accurate model of the weathered layer for good results. Frequently the weathered zones cannot be accurately modeled due to invalid initial assumptions about the zones"" properties or because the model is not well constrained by the first arrival time information. Also, traditional methods require the user to make accurate estimates of certain weathered layer attributes about which little is generally known. Traditional methods do not exploit the statistical independence of the offset term from the other terms employed in the present invention. Finally, traditional methods impose certain conditions or constraints on the model of the weathered layer which may not be consistent with the first arrival time information.
Analyzing the stacking responses of reflection data works better for estimating the short-period part of a statics correction term (so called residual statics), while estimation of the long-period statics correction term using reflection data alone is very unstable and inefficient. Therefore, analysis of refraction data for statics estimation correction solutions to reveal more information about the near surface, and correspondingly about statics, is customary in conventional seismic processing.
One well known prior art method of statics correction generation is the Generalized Reciprocal Method (GRM). The generalized reciprocal method has been widely applied to 2-D data (Palmer, D., 1981, An introduction to the generalized reciprocal method of seismic refraction interpretation: Geophysics, Soc. of Expl. Geophys., 46, 1508-1518). Unfortunately, for 3-D seismic the reciprocal method is difficult to apply due to the lack of reciprocal data. Another prior art method is called the Delay-Time method. The delay-time method assumes a near-surface model of locally flat layers on the scale of offset range. First arrivals were assumed to be the onset of head waves propagating along refracting interfaces of these layers. First arrival pick times are decomposed into delay times and refracting-layer velocities. Delay times are then converted to obtain layer thickness values and velocities assuming a critical angle of incidence on the refracting layers. An example of a delay time method is U.S. Pat. No. 6,424,920 to Osypov. The method uses differential delay times. Different types of refracted waves including head waves and diving waves are inverted. First, traveltime/offset functions for station in the seismic survey are estimated from first-arrival picks. Next, the traveltime/offset functions are transformed into velocity/depth functions to make up a near-surface model. Then, long-period statics are calculated using the derived near-surface model. Finally, short-period statics are estimated by a surface-consistent decomposition of the traveltime residuals.
The delay-time method is also used in other approaches, such as Generalized Linear Inversion (GLI) and head-wave refraction tomography (Hampson, D. and Russell, B., 1984, First-break interpretation using generalized linear inversion, 54th Ann. Internat. Mtg: Soc. of Expl. Geophys., Session:S10.4). In these methods, instead of the two-step inversion via delay times, traveltimes were inverted directly for layer thicknesses and velocities. However, these methods also had to deal with the problem of velocity/depth ambiguity. To address this issue, it was common practice to fix the weathering velocity prior to depth estimation. Incorporation of reflection data for joint inversion with refraction data could at times reduce the velocity/depth ambiguities. Unfortunately, near-surface reflections were very difficult to pick in normal production data.
One of the main problems with prior art refraction static correction methods has been that in areas with complex geology and rough terrains the simple model typically employed is not sufficient to explain important data features or adequately provide for good seismic imaging. Moreover, in order to fit the observed nonlinear move-out of first arrivals the offset range had to be limited or many more layers in the model had to be added. Both approaches made inversion even more unstable, ambiguous, time consuming and not sufficiently accurate.
Prior art solutions for static correction determinations are dependent upon input velocities and/or near-surface models that may be derived or estimated from weathering layers with highly variable velocities. This situation leads to inherent uncertainty whether the velocities are adequate for good statics correction solutions. It would be advantageous to have a method and system for statics correction solutions that did not depend upon user supplied input velocity parameters or models.
The present invention provides a method of determining source and receiver static corrections for seismic data that does not depend on user supplied velocity parameters or input models. The method comprises decomposing first break seismic data arrivals into source id terms, receiver id terms, offset terms and average datum correction terms. The source id terms and receiver id terms and offset terms are applied to the seismic data to obtain residual first break times. Average datum correction terms are determined by forming a nonlinear equation for a least-squares fit of the residual first break times as a function of elevation difference. At least one receiver static correction term is determined from the receiver id terms and receiver average datum correction terms. At least one source static correction is determined from said source id terms and source average datum correction terms.
The seismic data may be all within the same shot profile. The source, receiver offset terms and elevation difference terms may be determined using Gauss-Seidel iterations and discarding outlier values greater than a predetermined threshold. The nonlinear equation may be of the form y=Aexb, where y is residual time, x is the elevation difference from the average elevation, A is a scalar and b is exponent of the nonlinear best fit curve associated with the difference in elevation.